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ABSTRACT 

Aims. In order to understand the first stages of planet formation, when tiny grains aggregate to form planetesimals, one needs to 
simultaneously model grain growth, vertical settling and radial migration of dust in protoplanetary disks. In this study, we implement 
an analytical prescription for grain growth into a 3D two-phase hydrodynamics code to understand its effects on the dust distribution 
in disks. 

Methods. Following the analytic derivation of Stepinski & Valageas (1997), which assumes that grains stick perfectly upon collision, 
we implement a convenient and fast method of following grain growth in our 3D, two-phase (gas+dust) SPH code. We then follow 
the evolution of the size and spatial distribution of a dust population in a classical T Tauri star disk. 

Results. We find that the grains go through various stages of growth due to the complex interplay between gas drag, dust dynamics, and 
growth. Grains initially grow rapidly as they settle to the mid-plane, then experience a fast radial migration with little growth through 
the bulk of the disk, and finally pile-up in the inner disk where they grow more efficiently. This results in a bimodal distribution of 
grain sizes. Using this simple prescription of grain growth, we find that grains reach decimetric sizes in 10 5 years in the inner disk 
and survive the fast migration phase. 
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pH juttner & Yorkd [2Q0lt iDullemond & Domimkl Hjoj 
iTanaka et al] 120051: iNomura & Nakagawal 120061: ICieslal 120071) 
which describes the number density evolution of particles of a 
given mass range. The numerical solution of the Smoluchowski 
equation is challenging. Another approach is to use an analytic 
expression for the grain growth rate as a function of local d isk 
conditions (Stepinski & Valageas 1997; Haghighipour 2005|). In 
this study, we use this s econd approach and impleme nt the an- 
alytical prescription of Stepinski & Valagea| d 19971 hereafter 
SV97) in our 3D, two-phase (gas+dust) hydrodynamics code 
and follow the evolution of the grain size distribution in a proto- 
planetary disk. We validate our method on an axisymmetric disk 
here, before applying it to non-axisymmetric complex problems 
in future work. 



2. An analytical grain growth model 

|SW7j modeled the radial evolution of solid particles made of wa- 
ter ice in geometrically thin, turbulent, vertically isothermal pro- 
toplanetary disks. They described the gas and solid particle com- 
ponents as two separate phases coupled by aerodynamic forces 
(in the Epstein regime), and assumed that the evolution of the gas 
is unaffected by that of the solids. The particles can grow by co- 
agulation (they stick perfectly upon collision and therefore never 
shatter into smaller grains), evaporate or condense from vapor. 
Their size distribution at any radius and time is supposed to be 
narrowly peaked around a local mean value s(r, f). Self-gravity 
is neglected. 



1. Introduction 

The first steps of planet formation are governed by the build- 
up of planetesima ls due to the du st coagulation in proto- 
planetary disks (Dominik et al. 2007). Obser vational evidence 
for grain grow t h in d i sks is now common dApai et alj 12 004 ; 
Rodmann et al. 2006; Muzerolle et al 2006; iLommen et al.l 



l2007t lGraham et al.ll2.Q07D . Grains must grow from sub-yum sizes 
to planetesimal scale (kilometer size) objects in a fraction of the 
lifetime of the disk , which is estimated to be a few 10 7 years 
dHaisch et all 120011; ICarpenter etail 120051) . The timescales of 
grain growth, however, are unclear: some young disks show sig- 
natures of grain growth while old disks can show signatures of 
unprocessed grains and c oeval disks can show a rang e of grain 
sizes and dust processing ( Kessler-Silacci et al. 2006). 

Grains can grow via collisions and depending on their rel- 
ative velocity and on their che mical and physical properties 
dChokshi et al. 1993; Blum 2006), colliding grains can rebound, 
shatter or stick. Grains will settle vertica lly and migrate radially 
at different rates according to their size dWeiden schilling 19771 
iGaraud et al1l2004l iBarriere-Fouchet et alj|2005l) . leading to lo- 
cal density enhancements in the disk. Since grain growth is de- 
pendent on density, changes in the dust distribution will affect 
growth rates, which in turn will affect the d ynamics of the dust 
( Weidenschilling 1980; Haghighipour 2005). Therefore, growth, 
settling and migration need to be simulated together. 

Various models have been developed to describe the grain 
growth process. One approach is to use the time-dependent 
Smoluchowski coagulation equation (Weidenschilling 1980, 
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They obtained an analytic expression for the evolution of the 
particle size s given by 
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At 



Pi 



Sc 



(1) 



where Ro is the Rossby n umber for turbulent motions, a the 
Shak ura & Sunvaevl d 19731) viscosity parameter, pd the density 
of matter concentrated into solid particles, pd the bulk density 
of the grains, C s the local gas sound speed, and Sc the Schmidt 
number of the flow which estimates the effect of gas turbulence 
on the grains. They defined Sc as 



Sc = (i + n K Q 




(2) 



where £2k is the local keplerian velocity, f s the dust stopping 
time, v the mean relative vel ocity between gas and dus t, and V t a 
turbulent velocity. Note that lYoudin & Lithwickl d2007l) recently 
suggested that one should use Sc 1 + (Q#f ,) 2 instead. However 
it is not clear how this approximate value follows from their 
more complex exact expression (1 + (£2#fs) 2 ) 2 /(l +4(Q^f s ) 2 ). In 
this paper, we keep the formulation used by SV97 for the sake 
of self-consistency. 

ds 

The growth rate — depends on s via the stopping time 
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wh ere p g is the gas density. 

ISV97I were mostly interested in the radial distribution of 
solid particles and their sizes. However, in order to interpret 
observations of disks showing evidence of grain settling and 
growth, as well as to provide initial conditions for planet for- 
mation models, which are normally axisymmetric situations, one 
needs to know both the radial and vertical size and density distri- 
butions of grains. The complex interplay between the drag force, 
which causes solids to migrate radially and settle vertically, and 
the growth process makes a full numerical treatment of this prob- 
lem necessary. 



3. Grain growth in a 3D SPH code 

We have developed a 3D, two-phase (gas+dust) Smoothed 
Particles Hydrodynamics (SPH) code to model vertically 
isothermal, non self-gravitating protoplanetary disks. The two 
inter-penetrating phase s representing gas and dus t interact 
via aerodynamic drag. Barriere-Fouc het et al.l ([2005, hereafter 
BF05) describe the code and its limitations, and present the 
spatial distribution of dust grains ranging from 1 fim to 10 m 
in size resulting from radial migration and vertical settling. In 
thi s work , our aim is to implement the grain growth algorithm 
of ISV97I into our code and see how this simple prescription 
of grain growt h affe cts the dust dynamics by comparing with 
the results of BF05, While this could adequately be tested in 
2D, our ultimate goal is to study the observational signatures 
of grain growth in protoplanetary disks and apply the code to 
various non-axisymmetric problems, like disks with embedded 
planets. Thus we will be able to extend our previous work on 
the format ion of planetary gaps in the dust layers of p rotoplane- 
tary disks (Maddiso n et al.ll2007HFouchet et al.l l2007) to include 
grain growth, and study the stratifi cation of growing dust grains 
in disks (see, e.g JPinte et al.ll2007b . 



SPH is very well suited to the S V97 implementation of grain 
growth. As we do in our code, ISV97I describe gas and dust 
as two fluids, their disk configuration and thermodynamics are 
very similar to ours, and our SP H viscosity can be related to a 
Shakura-Sunyaev a viscosity dFouchet et al.l2007l) . Our test sim- 
ulations show that the gas disk is little affected by the evolution 
of the dust, as assu med b y lSV97l Therefore, it is straightforward 
to implement the ISV97I prescription of grain growth, given by 
Eq. (|T), in our code. 

Following the work of lSV97l all the dust particles must have 
the same initial size, so- The size s is then evolved using Eq. ([1) 
evaluated at the location of each SPH particle. We assume that 
the size s assigned to each SPH particle represents the typi- 
cal size of physical dust grains at its location in the dis k, at a 
given time. Again, this is very similar to the assumption of lSV97l 
for their local mean value s(r, t). Contrary to their work, we do 
not take evaporation into account: for our disk conditions (see 
Sect. 14. 11 1. ISV97I show that the evaporation radius is located be- 
tween 1 and 2 AU from the star and that after a million years only 
a small fraction of t he tota l solid material is evaporated due to 
migration. Similar to lSV97l our implementation of grain growth 
does not include the fragmentation of grains, even though it is 
likely to play an important role (see Sect. |6j. 

The mass of each SPH particle is kept constant to ensure 
kinetic energy and momentum conservation. As s can only in- 
crease, this implies that, over time, SPH dust particles repre- 
sent fewer but larger physical dust grains. However, even for 
km-sized planetesimals, each SPH particle still represents a very 
large number of physical particles, maintaining the validity of 
the numerical scheme. 



4. Simulations 

4.1. Setup 

We st udy grain growth in the typical T Tauri disk modeled in 
lBF05l with Mdisk = 0.02 M , composed of 99% gas and 1% dust 
by mass and orbiting a 1 M star. The dust grains have an in- 
trinsic density pd = 1 gem 3 . We choose an ini tial state for a 
gas di sk near equilibrium conditions. Following Havash i~et al.l 
(1985), we take the parameters of the Minimum Mass Solar 
Nebula for the radial dependence of the temperature (T oc r 3 ^ 4 ) 
and initial surface density (E oc r~ 3 ^ 2 ). The disk is locally isother- 
mal, i.e. the temperature follows a radial power law but is verti- 
cally constant. The sound speed then varies as C s oc r~ 3 ^ 8 and, 
given that H = C s /Qk, H/r varies as r 1 ^ . The disk is slightly 
flared with H/r = 0.05 at 100 AU. 

The smoothing length is computed by h\ oc (m i /p,) I/ ' 3 , where 
nij an d p,- are the mass and density of an SPH particle (see 
IBF05I) . The code will not be able to resolve a dust disk for which 
the semi-thickness is less than the smoothing length. However, 
turbulence in real disks stirs the solid particles and prevents 
the formation of such a thin dust layer. Indeed, the dust disk 
thickness reaches a steady state when th e settling and turbu- 



lent mixing are in equi libriu m (see, e.g.. iDubrulle et aljf l 995 ; 
iDullemond & DominiklfeOO^: ISchrapler & Henning|l20o4 . Our 



code does not yet include turbulent mixing and therefore can not 
reproduce this steady state. The SPH artificial viscosity terms 
are given by qtsph _ 0-1 and fepH = 0. 0, which ensures the 
corresponding Sha kura & Sunvaevl ( 1 1973b viscosity parameter 
a ~ .01 (as indicated by observations of p rotoplanetary disks 
— see lHartmann et al.lll998tlKing et al]|2007l) . 

We start with 200 000 gas particles that are distributed radi- 
ally so as to retrieve the expected power law for surface density 
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Fig. 1. Initial (top) and final (middle) dust density and final grain 
size distribution (bottom) in a meridian plane cut of the disk for 
So = lOyum. 



and randomly in the vertical direction because hydrostatic equi- 
librium is rapidly reached. The initial velocity of the gas parti- 
cles is Keplerian. Starting from this initial distribution, we allow 
the gas disk to relax for almost 8 000 years, which allows the 
pressure and artificial viscosity to smooth out the velocity field. 
Once the gas disk has relaxed, we then add an equal number of 
dust particles on top of the gas particles with the same velocity 
and allow the system to evolve. 

The disk extends initially from 20 to 300 AU. During the 
evolution, particles are removed from the simulation if they mi- 
grate inside of 20 AU and are assumed to be accreted by the 
star. The outer boundary is free and particles are only removed 
if they go beyond 400 AU. The viscous evolution of the disk is so 
slow that the gas surface density profile stays almost unchanged 
although the disk expands radially up to 400 AU. 

We ran a series of simulations with 400,000 SPH particles 
and initial grain sizes so ranging from 1 fim to 1 mm, in which 
the system is evolved for a total of 10 5 yr. The results presented 
in the next sections are time snapshots of the ongoing disk evo- 
lution. 

4.2. Results 

Figure Q] shows the density of the solid phase at dust injection 
(top) and at the end of the simulation (middle), along with the 
resulting size distribution (bottom) in a meridian plane cut of the 
disk for so = 10 yum. The vertical profiles of the gas and dust 
densities at the end of the simulation are shown in more detail 
in Fig. |2] The very efficient settling produces a thin dust disk, 
whereas the gas disk does not evolve on this timescale. 

Grain growth occurs very quickly: a few hundred years af- 
ter dust injection, the innermost grains have almost reached mm 
sizes (Fig. [3^), and a few thousand years later (Fig. [3j)), grain 
growth is visible over the entire disk, with sizes of a few cm in 
the inner region. This fast evolution leads to a radial size distri- 
bution showing a regular increase of grain size with decreasing 
distance from the star. At later times (Fig. [3J; and d), we see a 
change in the profile shape due to differences in migration effi- 
ciency in different parts of the disk. The overall distribution then 




-0.3 -0.2 -0.1 0.1 o.a 0.3 
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Fig. 2. Vertical (azimuthally averaged) profiles of the gas (thin 
lines) and dust (thick lines) densities at selected radii in the disk. 
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(a) t=8356 yr. 


(b) 1=15120 yr. 
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Fig. 3. Evolution of the radial grain size distribution for so = 
10 yi/m. 



evolves more slowly and shifts to larger sizes while keeping a 
roughly constant average slope. 

In Fig.|H we show the trajectories of seven individual parti- 
cles in the r — z plane for the so = 10 //m case. The top panel 
shows the particles settling to the mid-plane followed by their 
radial migration and the bottom panel shows the grain growth 
during this process. We see three stages: particles grow as they 
settle to the mid-plane and then start their radial migration (this 
behaviour, seen for all particles in Fig. [4] was also reported by 
Haghighipour (2005)), they then rapidly migrate inwards while 
growing very little (seen for P2, P3, P4, P5), and finally grow 
again while migrating very little (see for PI, P2, P3). The sec- 
ond stage of fast radial migration corresponds to the low slope 
seen in the central parts of the disk in Fig.[3j; and d: grains spend 
less time there and accumulate in the inner disk, where they grow 
more efficiently due to the higher density. 

The final size distribution (Fig. [3}l) shows a population of 
large grains over 10 cm in the inner disk, whereas the grains stay 
well below mm sizes in the outer disk. The histogram of final 
grain sizes for the so = 10 fim model (Fig.[5ji) shows a bimodal 
distribution. The first peak corresponds to the end of the first 
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Fig. 4. Trajectories of individual particles for so = 10 pm. 
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Fig. 5. Evolution of the histogram of grain sizes for so — 10 fxm. 



growth stage identified in Fig. [4] where grains reach a size of 
about 3 mm. The minimum around 8 mm is explained by the 
rapid migration of the second stage for grains of that size, which 
transports them to the very efficient growth region of the inner 
disk where they will populate the second peak around 10 cm. As 
is to be expected from Eq. ([T), the larger grains are found in the 
denser zones (see Fig.Q]). 

The histograms of grain sizes at different times displayed in 
Fig. [5] show that the signature of the rapid migration stage as a 
minimum around 8 mm appears only a few thousand years after 
the dust injection (Fig. [2)), as soon as the larger grains reach that 
size. It is visible thoughout the disk evolution, with a remarkably 
stable position (Fig. [5}: and d). 

Changing so from 1 pm to 1 mm only has an effect on the 
size distribution in the outer disk, where grains do not grow be- 
yond a few mm. In the rest of the disk, whatever their initial 
size, particles quickly reach the second stage and their subse- 
quent evolution is similar to the so = 10 pm case, leading to the 
same grain size distribution. 

In this paper, we have restricted our study to disks with a 
large inner radius of r- m = 20 AU. For lower values, simulations 
are much slower due to the smaller drag timestep in the central 
higher density regions. Simulations with r m = 3 AU were run 



for f = 20 000 yr and show the same size distribution outside of 
20 AU, with s continuing to increase as r decreases to 3 AU. 



5. Discussion 



Our results are in good agreement with those of lSV97l Indeed, 
with a disk of comparable mass to ours and an a viscosity pa- 
rameter of 10~ 2 , they found that grain growth is more efficient in 
the inner disk, where grains reach sizes of 1 cm at 20 AU after 
3 X 10 5 years. They explain their results by the combined action 
of both growth and migration: grains grow, reach sizes where 
they decouple from the gas and then migrate, thereby altering 
the dust surface density and grain size. Equilibrium occurs when 
the global flux of migrating particles at a given radius vanishes. 
However, they do not resolve the change of slope in the grain 
size distribution between 50 and 150 AU. Further, their inter- 
pretation remains qualitative and does not recover the different 
growth stages we observe. 

In order to interpret our results in a more quantitative way, 
we compare our simulations to an analytic expression of the re- 
sulting grain sizes derived under simplifying hypotheses. We 
first assume that v <k V t , i.e. that Sc is dominated by the ef- 
fect of gas-dust coupling, determined by the value of the non- 
dimensional stopping time T s = Qk t s . We introduce the non- 
dimensional grain size S = s/s opt , where 



^opt 



C sPg 
QltPd 



(4) 



is the optimal grai n size for radial migration (for details see 
iFouchet et al.l l2007). in order to write Sc =s 1 + S . Intere stingly, 
S — T s and the three drag regimes identified by BF05 are de- 
fined by S <K 1 (dust strongly coupled to the gas), S ~ 1 (most 
efficient dust settling and migration) and S » 1 (dust decoupled 
from the gas). 

To obtain an analytic solution for Eq. £[}, we then neglect the 
temporal variation of all quantities except s. Since ISV97I used 
constant values for Ro and a, the only remaining variables are 
Qk, p g , Pd and C s , which vary both in space and time. Assuming 
they are constant when following the evolution of the size of one 
particle amounts to assuming that, while it grows, the particle 
stays in the same position in the disk (therefore does not set- 
tle nor migrate) and that the disk structure does not evolve with 
time. While crude, this approximation is nonetheless useful in 
understanding the growth process and disentangling its effects 
from the complex interaction of those of the drag force. We de- 
fine the dimensionless time 



T = - +2 VSo 

T 



(■♦¥)■ 

to rewrite Eq. (fl} as 



with t = 



1 



Pg 



V2 3 /2RoaQ K Pd 



(5) 



dS_ 

dT 



l+S 
Its solution is 



(6) 



S = 



(8+9r 2 +3rVi6+9r 2 )' 



(8+9r 2 +3rVi6+9r 2 ) 



-2. (7) 
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t= 100268 yr , 









r (AU) 

Fig. 6. Radial grain size distribution computed from Eq. (0 at 
the end time of our simulations, for Sq - 10 /jm (black) and 
optimal size for radial migration s opt (grey). 
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showing a fast growth for small sizes and a slower growth for 
larger sizes, consistent with what was seen in Fig. [3] 

Although Eq. 01 associates one value of S to a given value 
of T, the set of normalization quantities (s opt , t) defining the non 
dimensional S and T is unique for each particle, and is computed 
with the disk parameters corresponding to its initial location. The 
resulting values of s will therefore be different for all particles. 

Starting from the same initial positions of dust particles of 
size sq, this simple model allows us, instead of evolving their 
size with the code, to directly compute with Eq. (0 the resulting 
size (i.e. the size computed for the end time of our simulations) 
of each particle in the approximation of a fixed position during 
growth. The resulting grain size distribution is shown in Fig. [6] 
for i() = 10 /urn. As s is proportional to p g via s opt , the observed 
spread in resulting grain size at a given r is caused by the spread 
in gas density for particles at different heights above the mid- 
plane. The resulting distribution reproduces the same grain size 
range as obtained at the end of the simulation, with decimetric 
grains in the inner disk, and has a slope close to that of the dis- 
tribution in Fig. [3]3. Of course, in that approximation, Eq. 01 
cannot reproduce the effect of migration and therefore does not 
show the changes in slope that it causes, which were seen in 
Fig. [3}; and d. 

The understanding of the growth process in the framework 
of our approximation which leads to Eq. (|6]l now allows us to 
reinterpret the results of grain growth in the presence of vertical 
settling and radial migration. For the range of sq values used in 
our simulations, so < s opl over most of the disk (see Fig. [6), and 
the grains are initially in the first regime of fast growth (S <K 
1) identified in Sect. 14.21 and found in Eq. {3), corresponding 
to the strongly coupled drag regime, with very little migration. 
Outside of ~ 150 AU, resulting grain sizes stay below s opt in the 
simulated timescale and grains do not reach the rapid migration 
stage. Indeed, the growth timescale r oc Q" 1 oc r k is longer in 
the outer disk. Inside of this radius, sizes do reach s opt (i.e. S ~ 1, 
corresponding to the intermediate drag regime) and grains go 
through a phase of fast inwards migration, explaining the low 
slope in the central region of Fig. 01 and the stable location of the 
minimum in the histograms of grain sizes (Fig. [5J5— d). Indeed, 



i opt depends on variables describing the gas disk, which does 
not evolve in the timescale of our simulations. Once they have 
grown to sizes greater than s opt , the migration efficiency drops 
as dust starts to decouple from the gas (in the weak drag regime, 
defined by S » 1) and grains pile up inside of ~ 50 AU and 
continue to grow there, more slowly. 

This result differs from that of lWeidenschillinj d 1977b . who 
models the particular case of the Minimum Mass Solar Nebula. 
With his disk parameters, he found that the grains with the 
largest radial migration velocity were meter-sized and that their 
consequent survival time in the nebula was much s horter than 
the di sk lifetime. This led to the persisting (see, e.g jNatta et al.l 
2007) idea that growing grains which reached meter sizes rapidly 
fall onto the star (or more likely are evaporated in the inner disk), 
thus resulting in a potential problem for the planet formation 
process. This so-called "meter-size barrier" could be overcome 
if grains could grow to larger sizes within only a few hundred 
years, which would be very difficult accord ing to the c urrent 
understanding of solid particle aggregation dBluml 120061) . Our 
simulations show that, in our slightly more massive nebula rep- 
resentative of many observed T Tauri disks, growing grains can 
survive the fast migration stage (occurring for ~ 8 mm grains in 
our disk conditions) and reach larger sizes with longer survival 
times, allowing for the possibility of planet formation inside of 
-50AU. 



6. Conclusion 

We have implemented a mechanism able to treat grain growth in 
protoplanetary disks via the analytical expression of lSV97l into 
our two-phase SPH code. We simulated for the first time the full 
3D evolution of a typical T Tauri disk, following the simultane- 
ous radial migration, vertical settling, and growth of solid parti- 
cles. Their interplay is complex: dynamics affects grain growth 
by modifying local physical quantities such as density or relative 
velocity. Conversely, grain growth also acts on dust dynamics: 
where non-growing grains would either stay well mixe d wit h 
the gas or settle and migrate according to their sizes ([BF05), 
growing grains will go through various stages and produce to- 
tally different spatial distributions. They initially grow rapidly 
as they settle to the mid-plane, then experience a fast radial mi- 
gration with little growth through the bulk of the disk, and finally 
pile-up in the inner disk where they grow more efficiently. This 
results in a bimodal distribution of grain sizes, with the largest 
grains found in the denser inner disk, where growth is most ef- 
ficient. The survival times of the solid particles are longer than 
previously found, which has an implication on planet formation. 

We find that grains grow very quickly: they reach decimet- 
ric size in 10 5 yr. This is i n gene ral agreement with the results 
of iDullemond & Dominikl d2005l) where shattering is neglected. 
They used the different approach of solving the Smoluchowski 
equation to study the growth of settling, but non-migrating, dust 
and in 10 5 yr formed grains of maximum sizes ranging from 
1 cm to over 1 m depending on their model parameters. Similar 
to them, we also find that the small grains only survive in the 
very outer disk and are depleted too rapidly elsewhere to be con- 
sistent with infrared observations of disks, highlighting the im- 
portance of shattering. 

In order to compute synthetic images from our simulations 
and compare them to the observations, one would have to as- 
sume that the collisional cascade resulting from the inclusion of 
shattering would produce a whole particle distribution from the 
maximum size at a given radius shown in Fig. [3jl down to sub- 
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jum size, desc ribed by a quasi-steady power law as argued by 
lGaraudld2007h . 

The method we used to treat grain growth can easily be ap- 
plied to other analytical prescriptions. The development of a 
more detailed model is necessary for a realistic description of 
grain growth in protoplanetary disks. In addition to shattering, 
one needs to take into account other processes such as micro- 
scopic interactions between the grains, kinetic energy dissipation 
and grain porosity. This is the subject of a forthcoming paper. 
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